##############################################################
##############################################################
####### Replication code for "Do museums promote reconciliation? Evidence from a field experiment," Journal of Politics
####### This file includes code for all analyses and figures in the online appendix
##############################################################
##############################################################

require("sandwich")
require("plyr")
require("lmtest")
require(dplyr)
require(gridExtra)
require("RColorBrewer")
require(ggplot2)

##############################################################
###### Read in data and establish main functions
##############################################################
load(file = "all.Rdata")

### ATE FUNCTIONS ##
# This estimates ATE when we have a pre-treatment measurement
est.ate<-function(dv, predv, df){
  predv <- f(predv)
  dv <- f(dv)
  summary(fit.1 <- lm(dv~df$treat + df$pre_ideology_1 +
                        + predv*df$date_diff + df$base_gender +df$age + df$v))
  vcv <- vcovHC(fit.1)
  n <- nobs(fit.1)
  result <- coeftest(fit.1, vcv)[2, 1:4] / itt.d
  result <- list("point.estimate" = result[1],"se"=result[2],"pvalue"=result[4], "obs" = n)
  return(result)
}
# This estimates ATE when we don't have a pre-treatment measurement
est.ate.np<-function(dv, df){
  dv <- f(dv)
  summary(fit.1 <- lm(dv~df$treat + df$pre_ideology_1 + df$base_gender +df$age+df$v))
  vcv <- vcovHC(fit.1)
  n <- nobs(fit.1)
  result <- coeftest(fit.1, vcv)[2, 1:4] / itt.d
  result <- list("point.estimate" = result[1],"se"=result[2],"pvalue"=result[4], "obs" = n)
  return(result)
}

# this function recodes NAs to the mean, per our PAP
f <- function(x){ 
  m <- mean(x, na.rm = TRUE) 
  x[is.na(x)] <- m 
  x 
}
## recode covariates to means
all$age <- f(all$age)
all$pre_ideology_1 <- f(all$pre_ideology_1)
all$base_gender <- f(all$base_gender)
all$date_diff <- f(all$date_diff)
# split dataset into left, right, and related to victim for heterogeneous analyses
left <- all[all$right == 0,]
right <- all[all$right == 1,]
itt.d <- all$itt.d

##############################################################
###### Table A1: Number of respondents by condition
##############################################################

# Note: we do not have data on those who did not show up or opted not to comply with their random assignment

# completed by condition
sum(all$treat==1)
sum(all$treat==0)

# first follow up by condition
sum(all$f1[all$treat==1])
sum(all$f1[all$treat==0])

# second follow up by condition
sum(all$f2[all$treat==1])
sum(all$f2[all$treat==0])

# third follow up by condition
sum(all$f3[all$treat==1])
sum(all$f3[all$treat==0])

##############################################################
###### Table A2: Covariate balance
##############################################################
t.test(all$age~all$treat)
t.test(all$base_gender~all$treat)
t.test(all$pre_ideology_1~all$treat)
t.test(all$v~all$treat)
t.test(all$pre_political_interest~all$treat)
t.test(all$pre_party_id~all$treat)
t.test(all$pre_positive~all$treat)
t.test(all$pre_negative~all$treat)
t.test(all$pre_conf_gov~all$treat)

##############################################################
###### Table A3: Perceptions of museum by ideology
##############################################################
t.test(all$mm_obj~all$right)
t.test(all$mm_views_like~all$right)
t.test(all$mm_views_content~all$right)
t.test(all$mm_views_inhibit~all$right)
t.test(all$mm_views_important~all$right)
t.test(all$mm_new~all$right)

##############################################################
###### Table A5: Perceptions of inequality after visiting the MMDH
##############################################################
est.ate(all$current_ineq, all$pre_current_ineq, all)

##############################################################
###### Table A6: Full regression results, Political institutions
##############################################################

# See lines 64-72 in "BPV_museums_maintext.R"

##############################################################
###### Table A7: Full regression results, Political institutions by ideology
##############################################################

# See lines 75-93 in "BPV_museums_maintext.R"
# Interactions (final three columns of table) reproduced here
est.ate.int<-function(dv, predv, df){
  predv <- f(predv)
  dv <- f(dv)
  summary(fit.1 <- lm(dv~df$treat*df$pre_ideology_1+
                        + predv*df$date_diff + df$base_gender +df$age + df$v))
  vcv <- vcovHC(fit.1)
  n <- nobs(fit.1)
  result <- coeftest(fit.1, vcv)[9, 1:4] / itt.d
  result <- list("point.estimate" = result[1],"se"=result[2],"pvalue"=result[4], "obs" = n)
  return(result)
}

est.ate.int(all$democracy, all$pre_democracy, all)
est.ate.int(all$military_gov, all$pre_military_gov, all)
est.ate.int(all$inst_gov, all$pre_inst_gov, all)
est.ate.int(all$inst_mil, all$pre_inst_mil, all)
est.ate.int(all$inst_police, all$pre_inst_police, all)
est.ate.int(all$conf_gov, all$pre_conf_gov, all)
est.ate.int(all$conf_mil, all$pre_conf_mil, all)
est.ate.int(all$conf_police, all$pre_conf_police, all)
est.ate.int(all$conf_church, all$pre_conf_church, all)

##############################################################
###### Table A8: Full regression results, Transitional justice
##############################################################

# See lines 164-171 in "BPV_museums_maintext.R"

##############################################################
###### Table A9: Full regression results, Transitional justice by ideology
##############################################################

# See lines 75-93 in "BPV_museums_maintext.R"
# Interactions (final three columns of table) reproduced below
est.ate.int.np<-function(dv, df){
  dv <- f(dv)
  summary(fit.1 <- lm(dv~df$treat*df$pre_ideology_1+
                        + df$base_gender +df$age + df$v))
  vcv <- vcovHC(fit.1)
  n <- nobs(fit.1)
  result <- coeftest(fit.1, vcv)[7, 1:4] / itt.d
  result <- list("point.estimate" = result[1],"se"=result[2],"pvalue"=result[4], "obs" = n)
  return(result)
}

est.ate.int.np(all$justice_advance, all)
est.ate.int.np(all$justice_account, all)
est.ate.int(all$current_recomp, all$pre_current_recomp, all)
est.ate.int.np(all$tj_judicial, all)
est.ate.int.np(all$tj_apology, all)
est.ate.int.np(all$policies_apologize, all)
est.ate.int.np(all$policies_compensate, all)
est.ate.int.np(all$policies_pardon, all)

##############################################################
###### Table A10: Full regression results, emotions
##############################################################

# See lines 261-297 in "BPV_museums_maintext.R"

##############################################################
###### Table A11: Number of respondents by condition
##############################################################

# See code for A1, above

##############################################################
###### Table A12: Test for differential attrition
##############################################################

## F-test
diff.att.full <- lm(all$observerd~(all$treat*all$right)
                    +(all$treat*all$base_gender) + (all$treat*all$age) + 
                      (all$treat*all$v))
diff.att.1.full <- lm(all$f1~(all$treat*all$right)
                      +(all$treat*all$base_gender) + (all$treat*all$age) +
                        (all$treat*all$v))
diff.att.2.full <- lm(all$f2~(all$treat*all$right)
                      +(all$treat*all$base_gender) + (all$treat*all$age) +
                        (all$treat*all$v))
diff.att.3.full <- lm(all$f3~(all$treat*all$right)
                      +(all$treat*all$base_gender) + (all$treat*all$age) +
                        (all$treat*all$v))
diff.att.red <- lm(all$observerd~all$treat+all$right
                   +all$base_gender+all$age + all$v)
diff.att.1.red <- lm(all$f1~all$treat+all$right
                     +all$base_gender+all$age + all$v)
diff.att.2.red <- lm(all$f2~all$treat+all$right
                     +all$base_gender+all$age + all$v)
diff.att.3.red <- lm(all$f3~all$treat+all$right
                     +all$base_gender+all$age + all$v)
anova(diff.att.full, diff.att.red)
anova(diff.att.1.full, diff.att.1.red)
anova(diff.att.2.full, diff.att.2.red)
anova(diff.att.3.full, diff.att.3.red)

##############################################################
###### Table A13: Differential attrition by pre-treatment covariates
##############################################################

#age
summary(lm(f(all$f1)~f(all$age)))
summary(lm(f(all$f2)~f(all$age)))
summary(lm(f(all$f3)~f(all$age)))

# gender
summary(lm(f(all$f1)~f(all$base_gender)))
summary(lm(f(all$f2)~f(all$base_gender)))
summary(lm(f(all$f3)~f(all$base_gender)))

# ideology
summary(lm(f(all$f1)~f(all$pre_ideology_1)))
summary(lm(f(all$f2)~f(all$pre_ideology_1)))
summary(lm(f(all$f3)~f(all$pre_ideology_1)))

# economic situation
summary(lm(f(all$f1)~f(all$pre_economic_situation)))
summary(lm(f(all$f2)~f(all$pre_economic_situation)))
summary(lm(f(all$f3)~f(all$pre_economic_situation)))

# political interest
summary(lm(f(all$f1)~f(all$pre_political_interest)))
summary(lm(f(all$f2)~f(all$pre_political_interest)))
summary(lm(f(all$f3)~f(all$pre_political_interest)))

# religiosity
summary(lm(f(all$f1)~f(all$pre_religion_importance)))
summary(lm(f(all$f2)~f(all$pre_religion_importance)))
summary(lm(f(all$f3)~f(all$pre_religion_importance)))

# museum visits
summary(lm(f(all$f1)~f(all$totalmuseums)))
summary(lm(f(all$f2)~f(all$totalmuseums)))
summary(lm(f(all$f3)~f(all$totalmuseums)))

# trust in the government
summary(lm(f(all$f1)~f(all$pre_conf_gov)))
summary(lm(f(all$f2)~f(all$pre_conf_gov)))
summary(lm(f(all$f3)~f(all$pre_conf_gov)))

# satisfaction with the government
summary(lm(f(all$f1)~f(all$pre_inst_gov)))
summary(lm(f(all$f2)~f(all$pre_inst_gov)))
summary(lm(f(all$f3)~f(all$pre_inst_gov)))

# inequality is a problem
summary(lm(f(all$f1)~f(all$pre_current_ineq)))
summary(lm(f(all$f2)~f(all$pre_current_ineq)))
summary(lm(f(all$f3)~f(all$pre_current_ineq)))

# positive emotions
summary(lm(f(all$f1)~f(all$pre_positive)))
summary(lm(f(all$f2)~f(all$pre_positive)))
summary(lm(f(all$f3)~f(all$pre_positive)))

#negative emotions
summary(lm(f(all$f1)~f(all$pre_negative)))
summary(lm(f(all$f2)~f(all$pre_negative)))
summary(lm(f(all$f3)~f(all$pre_negative)))

##############################################################
###### Table A14: Differential attrition by round 1 responses
##############################################################

# pol institutions index
summary(lm(f(all$f1)~f(all$pol.inst.index)))
summary(lm(f(all$f2)~f(all$pol.inst.index)))
summary(lm(f(all$f3)~f(all$pol.inst.index)))

# military gov
summary(lm(f(all$f1)~f(all$military_gov)))
summary(lm(f(all$f2)~f(all$military_gov)))
summary(lm(f(all$f3)~f(all$military_gov)))

# tj index
summary(lm(f(all$f1)~f(all$tj.index)))
summary(lm(f(all$f2)~f(all$tj.index)))
summary(lm(f(all$f3)~f(all$tj.index)))

# compensation
summary(lm(f(all$f1)~f(all$current_recomp)))
summary(lm(f(all$f2)~f(all$current_recomp)))
summary(lm(f(all$f3)~f(all$current_recomp)))

# pardon
summary(lm(f(all$f1)~f(all$policies_pardon)))
summary(lm(f(all$f2)~f(all$policies_pardon)))
summary(lm(f(all$f3)~f(all$policies_pardon)))

# negative emotions
summary(lm(f(all$f1)~f(all$negative)))
summary(lm(f(all$f2)~f(all$negative)))
summary(lm(f(all$f3)~f(all$negative)))

# positive emotions
summary(lm(f(all$f1)~f(all$positive)))
summary(lm(f(all$f2)~f(all$positive)))
summary(lm(f(all$f3)~f(all$positive)))


##############################################################
###### Table A15: Political institutions adjusted for multiple comparisons
##############################################################

# conducted using results from "BPV_museums_maintext.R" and EGAP calculator - https://egap.shinyapps.io/multiple-comparisons-app/
# adjusted significance with Bejamini and Hochberg correction

##############################################################
###### Table A16: Transitional justice adjusted for multiple comparisons
##############################################################

# conducted using results from "BPV_museums_maintext.R" and EGAP calculator - https://egap.shinyapps.io/multiple-comparisons-app/
# adjusted significance with Bejamini and Hochberg correction

##############################################################
######Table A17. General museum impressions by recoded ideology.
##############################################################

# Split up by ideology
## RECODE FOR ROBUSTNESS HERE ##
all$pre_ideology_1 <- f(all$pre_ideology_1)
all$right <- ifelse(all$pre_ideology_1 > 5, 1,0)
left <- all[all$right == 0,]
right <- all[all$right == 1,]

# mean values on dvs
t.test(all$mm_obj~all$right)
t.test(all$mm_views_like~all$right)
t.test(all$mm_views_content~all$right)
t.test(all$mm_views_inhibit~all$right)
t.test(all$mm_views_important~all$right)
t.test(all$mm_new~all$right)

##############################################################
######Table A18. Political institutions by recoded ideology
##############################################################

dem.right <- est.ate(right$democracy, right$pre_democracy, right)
mil.right <- est.ate(right$military_gov, right$pre_military_gov, right)
gov_sat.right <- est.ate(right$inst_gov, right$pre_inst_gov, right)
mil_sat.right <- est.ate(right$inst_mil, right$pre_inst_mil, right)
pol_sat.right <- est.ate(right$inst_police, right$pre_inst_police, right)
gov_trust.right <- est.ate(right$conf_gov, right$pre_conf_gov, right)
mil_trust.right <- est.ate(right$conf_mil, right$pre_conf_mil, right)
pol_trust.right <- est.ate(right$conf_police, right$pre_conf_police, right)
church_trust.right <- est.ate(right$conf_church, right$pre_conf_church, right)

dem.left <- est.ate(left$democracy, left$pre_democracy, left)
mil.left <- est.ate(left$military_gov, left$pre_military_gov, left)
gov_sat.left <- est.ate(left$inst_gov, left$pre_inst_gov, left)
mil_sat.left <- est.ate(left$inst_mil, left$pre_inst_mil, left)
pol_sat.left <- est.ate(left$inst_police, left$pre_inst_police, left)
gov_trust.left <- est.ate(left$conf_gov, left$pre_conf_gov, left)
mil_trust.left <- est.ate(left$conf_mil, left$pre_conf_mil, left)
pol_trust.left <- est.ate(left$conf_police, left$pre_conf_police, left)
church_trust.left <- est.ate(left$conf_church, left$pre_conf_church, left)
## interactions (for appendix)

dem_int <- est.ate.int(all$democracy, all$pre_democracy, all)
mil.int <- est.ate.int(all$military_gov, all$pre_military_gov, all)
gov_sat.int <- est.ate.int(all$inst_gov, all$pre_inst_gov, all)
mil_sat.int <- est.ate.int(all$inst_mil, all$pre_inst_mil, all)
pol_sat.int <- est.ate.int(all$inst_police, all$pre_inst_police, all)
gov_trust.int <- est.ate.int(all$conf_gov, all$pre_conf_gov, all)
mil_trust.int <- est.ate.int(all$conf_mil, all$pre_conf_mil, all)
pol_trust.int <- est.ate.int(all$conf_police, all$pre_conf_police, all)
church_trust.int <- est.ate.int(all$conf_church, all$pre_conf_church, all)

##############################################################
######Table A19. Transitional justice by recoded ideology
##############################################################

advance.right <- est.ate.np(right$justice_advance, right)
justice_account.right <- est.ate.np(right$justice_account, right)
compensation.right <- est.ate(right$current_recomp, right$pre_current_recomp, right)
judicial.right <- est.ate.np(right$tj_judicial, right)
inst_apology.right <- est.ate.np(right$tj_apology, right)
apologize.right <- est.ate.np(right$policies_apologize, right)
compensate.right <- est.ate.np(right$policies_compensate, right)
pardoned.right <- est.ate.np(right$policies_pardon, right)

advance.left <- est.ate.np(left$justice_advance, left)
justice_account.left <- est.ate.np(left$justice_account, left)
compensation.left <- est.ate(left$current_recomp, left$pre_current_recomp, left)
judicial.left <- est.ate.np(left$tj_judicial, left)
inst_apology.left <- est.ate.np(left$tj_apology, left)
apologize.left <- est.ate.np(left$policies_apologize, left)
compensate.left <- est.ate.np(left$policies_compensate, left)
pardoned.left <- est.ate.np(left$policies_pardon, left)

##############################################################
###### Figure A3: Persistence of responses across treatment groups and ideologies
##############################################################

# unrecode ideology
# Split up by ideology
## RECODE FOR ROBUSTNESS HERE ##
all$pre_ideology_1 <- f(all$pre_ideology_1)
all$right <- ifelse(all$pre_ideology_1 > 4, 1,0)
left <- all[all$right == 0,]
right <- all[all$right == 1,]

est.ate.np.f<-function(dv){
  summary(fit.1 <- lm(dv~all$treat + all$pre_ideology_1 + all$base_gender +all$age+all$v))
  vcv <- vcovHC(fit.1)
  n <- nobs(fit.1)
  result <- coeftest(fit.1, vcv)[2, 1:4] / itt.d
  result <- list("point.estimate" = result[1],"se"=result[2],"pvalue"=result[4], "obs" = n)
  return(result)
}

## PARDON
pardoned <- est.ate.np(all$policies_pardon, all)
pardoned_f1<- est.ate.np.f(all$policies_pardon_f1)
pardoned_f2 <- est.ate.np.f(all$policies_pardon_f2)
pardoned_f3 <- est.ate.np.f(all$policies_pardon_f3)
prop.table(table(all$policies_pardon[all$treat==1&all$right==1]))
prop.table(table(all$policies_pardon_f1[all$treat==1&all$right==1]))
prop.table(table(all$policies_pardon_f2[all$treat==1&all$right==1]))
prop.table(table(all$policies_pardon_f3[all$treat==1&all$right==1]))

all$pardon <- ifelse(all$policies_pardon==1|all$policies_pardon==0,0,1)
all$pardon_f1 <- ifelse(all$policies_pardon_f1==1|all$policies_pardon_f1==0,0,1)
all$pardon_f2 <- ifelse(all$policies_pardon_f2==1|all$policies_pardon_f2==0,0,1)
all$pardon_f3 <- ifelse(all$policies_pardon_f3==1|all$policies_pardon_f3==0,0,1)

pardon.df <- data.frame(all$pardon, all$pardon_f1,all$pardon_f2,all$pardon_f3,all$treat,all$right)
pardon.df.treat.right <- pardon.df[pardon.df$all.treat==1&pardon.df$all.right==1,]
pardon.df.treat.right$num <- rowSums(pardon.df.treat.right==1)-2  
prop.table(table(pardon.df.treat.right$num))

pardon.df.treat.left <- pardon.df[pardon.df$all.treat==1&pardon.df$all.right==0,]
pardon.df.treat.left$num <- rowSums(pardon.df.treat.left==1)-1  
prop.table(table(pardon.df.treat.left$num))

pardon.df.control.right <- pardon.df[pardon.df$all.treat==0&pardon.df$all.right==1,]
pardon.df.control.right$num <- rowSums(pardon.df.control.right==1)-1  
prop.table(table(pardon.df.control.right$num))

pardon.df.control.left <- pardon.df[pardon.df$all.treat==0&pardon.df$all.right==0,]
pardon.df.control.left$num <- rowSums(pardon.df.control.left==1)  
pardon.df.control.left$num <- factor(pardon.df.control.left$num, levels = c(0:4))
prop.table(table(pardon.df.control.left$num))
df <- data.frame(waveschosen=c(0:4),percent=NA)

df$percent <- prop.table(table(pardon.df.control.left$num))
df$group <- "control left"
df2 <- data.frame(waveschosen=c(0:4),percent=NA)
pardon.df.control.right$num <- factor(pardon.df.control.right$num, levels = c(0:4))

df2$percent <- prop.table(table(pardon.df.control.right$num))
df2$group <- "control right"

df3 <- data.frame(waveschosen=c(0:4),percent=NA)
pardon.df.treat.right$num <- factor(pardon.df.treat.right$num, levels = c(0:4))

df3$percent <- prop.table(table(pardon.df.treat.right$num))
df3$group <- "treat right"

df4 <- data.frame(waveschosen=c(0:4),percent=NA)
pardon.df.treat.left$num <- factor(pardon.df.treat.left$num, levels = c(0:4))

df4$percent <- prop.table(table(pardon.df.treat.left$num))
df4$group <- "treat left"

pardon.df <- rbind (df,df2,df3,df4)

## trust church
all$pre_conf_church <- ifelse(all$pre_conf_church==1|all$pre_conf_church==0,0,1)
all$conf_church <- ifelse(all$conf_church==1|all$conf_church==0,0,1)
all$conf_church_f1 <- ifelse(all$conf_church_f1==1|all$conf_church_f1==0,0,1)
all$conf_church_f2 <- ifelse(all$conf_church_f2==1|all$conf_church_f2==0,0,1)
all$conf_church_f3 <- ifelse(all$conf_church_f3==1|all$conf_church_f3==0,0,1)


trust_church.df <- data.frame(all$conf_church, all$conf_church_f1,all$conf_church_f2,all$conf_church_f3,all$treat,all$right)
trust_church.treat.right <- trust_church.df[trust_church.df$all.treat==1&trust_church.df$all.right==1,]
trust_church.treat.right$num <- rowSums(trust_church.treat.right==1)-2  
prop.table(table(trust_church.treat.right$num))

trust_church.treat.left <- trust_church.df[trust_church.df$all.treat==1&trust_church.df$all.right==0,]
trust_church.treat.left$num <- rowSums(trust_church.treat.left==1)-1  
prop.table(table(trust_church.treat.left$num))

trust_church.control.right <- trust_church.df[trust_church.df$all.treat==0&trust_church.df$all.right==1,]
trust_church.control.right$num <- rowSums(trust_church.control.right==1)-1  
prop.table(table(trust_church.control.right$num))

trust_church.control.left <- trust_church.df[trust_church.df$all.treat==0&trust_church.df$all.right==0,]
trust_church.control.left$num <- rowSums(trust_church.control.left==1)  
trust_church.control.left$num <- factor(trust_church.control.left$num, levels = c(0:4))
prop.table(table(trust_church.control.left$num))
df <- data.frame(waveschosen=c(0:4),percent=NA)

df$percent <- prop.table(table(trust_church.control.left$num))
df$group <- "control left"
df2 <- data.frame(waveschosen=c(0:4),percent=NA)
trust_church.control.right$num <- factor(trust_church.control.right$num, levels = c(0:4))

df2$percent <- prop.table(table(trust_church.control.right$num))
df2$group <- "control right"

df3 <- data.frame(waveschosen=c(0:4),percent=NA)
trust_church.treat.right$num <- factor(trust_church.treat.right$num, levels = c(0:4))

df3$percent <- prop.table(table(trust_church.treat.right$num))
df3$group <- "treat right"

df4 <- data.frame(waveschosen=c(0:4),percent=NA)
trust_church.treat.left$num <- factor(trust_church.treat.left$num, levels = c(0:4))

df4$percent <- prop.table(table(trust_church.treat.left$num))
df4$group <- "treat left"

churchtrust.df <- rbind (df,df2,df3,df4)

## satisfaction with government
all$pre_inst_gov <- ifelse(all$pre_inst_gov==1|all$pre_inst_gov==0,0,1)
all$inst_gov <- ifelse(all$inst_gov==1|all$inst_gov==0,0,1)
all$inst_gov_f1 <- ifelse(all$inst_gov_f1==1|all$inst_gov_f1==0,0,1)
all$inst_gov_f2 <- ifelse(all$inst_gov_f2==1|all$inst_gov_f2==0,0,1)
all$inst_gov_f3 <- ifelse(all$inst_gov_f3==1|all$inst_gov_f3==0,0,1)

inst_gov.df <- data.frame(all$inst_gov, all$inst_gov_f1,all$inst_gov_f2,all$inst_gov_f3,all$treat,all$right)
inst_gov.treat.right <- inst_gov.df[inst_gov.df$all.treat==1&inst_gov.df$all.right==1,]
inst_gov.treat.right$num <- rowSums(inst_gov.treat.right==1)-2  
prop.table(table(inst_gov.treat.right$num))

inst_gov.treat.left <- inst_gov.df[inst_gov.df$all.treat==1&inst_gov.df$all.right==0,]
inst_gov.treat.left$num <- rowSums(inst_gov.treat.left==1)-1  
prop.table(table(inst_gov.treat.left$num))

inst_gov.control.right <- inst_gov.df[inst_gov.df$all.treat==0&inst_gov.df$all.right==1,]
inst_gov.control.right$num <- rowSums(inst_gov.control.right==1)-1  
prop.table(table(inst_gov.control.right$num))

inst_gov.control.left <- inst_gov.df[inst_gov.df$all.treat==0&inst_gov.df$all.right==0,]
inst_gov.control.left$num <- rowSums(inst_gov.control.left==1)  
inst_gov.control.left$num <- factor(inst_gov.control.left$num, levels = c(0:4))
prop.table(table(inst_gov.control.left$num))
df <- data.frame(waveschosen=c(0:4),percent=NA)

df$percent <- prop.table(table(inst_gov.control.left$num))
df$group <- "control left"
df2 <- data.frame(waveschosen=c(0:4),percent=NA)
inst_gov.control.right$num <- factor(inst_gov.control.right$num, levels = c(0:4))

df2$percent <- prop.table(table(inst_gov.control.right$num))
df2$group <- "control right"

df3 <- data.frame(waveschosen=c(0:4),percent=NA)
inst_gov.treat.right$num <- factor(inst_gov.treat.right$num, levels = c(0:4))

df3$percent <- prop.table(table(inst_gov.treat.right$num))
df3$group <- "treat right"

df4 <- data.frame(waveschosen=c(0:4),percent=NA)
inst_gov.treat.left$num <- factor(inst_gov.treat.left$num, levels = c(0:4))

df4$percent <- prop.table(table(inst_gov.treat.left$num))
df4$group <- "treat left"

govsat.df <- rbind (df,df2,df3,df4)



## satisfaction with democracy

all$pre_democracy <- ifelse(all$pre_democracy==1|all$pre_democracy==0,0,1)
all$democracy <- ifelse(all$democracy==1|all$democracy==0,0,1)
all$democracy_f1 <- ifelse(all$democracy_f1==1|all$democracy_f1==0,0,1)
all$democracy_f2 <- ifelse(all$democracy_f2==1|all$democracy_f2==0,0,1)
all$democracy_f3 <- ifelse(all$democracy_f3==1|all$democracy_f3==0,0,1)


democracy.df <- data.frame(all$democracy, all$democracy_f1,all$democracy_f2,all$democracy_f3,all$treat,all$right)
democracy.treat.right <- democracy.df[democracy.df$all.treat==1&democracy.df$all.right==1,]
democracy.treat.right$num <- rowSums(democracy.treat.right==1)-2  
prop.table(table(democracy.treat.right$num))

democracy.treat.left <- democracy.df[democracy.df$all.treat==1&democracy.df$all.right==0,]
democracy.treat.left$num <- rowSums(democracy.treat.left==1)-1  
prop.table(table(democracy.treat.left$num))

democracy.control.right <- democracy.df[democracy.df$all.treat==0&democracy.df$all.right==1,]
democracy.control.right$num <- rowSums(democracy.control.right==1)-1  
prop.table(table(democracy.control.right$num))

democracy.control.left <- democracy.df[democracy.df$all.treat==0&democracy.df$all.right==0,]
democracy.control.left$num <- rowSums(democracy.control.left==1)  
democracy.control.left$num <- factor(democracy.control.left$num, levels = c(0:4))
prop.table(table(democracy.control.left$num))
df <- data.frame(waveschosen=c(0:4),percent=NA)

df$percent <- prop.table(table(democracy.control.left$num))
df$group <- "control left"
df2 <- data.frame(waveschosen=c(0:4),percent=NA)
democracy.control.right$num <- factor(democracy.control.right$num, levels = c(0:4))

df2$percent <- prop.table(table(democracy.control.right$num))
df2$group <- "control right"

df3 <- data.frame(waveschosen=c(0:4),percent=NA)
democracy.treat.right$num <- factor(democracy.treat.right$num, levels = c(0:4))

df3$percent <- prop.table(table(democracy.treat.right$num))
df3$group <- "treat right"

df4 <- data.frame(waveschosen=c(0:4),percent=NA)
democracy.treat.left$num <- factor(democracy.treat.left$num, levels = c(0:4))

df4$percent <- prop.table(table(democracy.treat.left$num))
df4$group <- "treat left"

democracy.df <- rbind (df,df2,df3,df4)

## military government

military_gov.df <- data.frame(all$military_gov, all$military_gov_f1,all$military_gov_f2,all$military_gov_f3,all$treat,all$right)
military_gov.treat.right <- military_gov.df[military_gov.df$all.treat==1&military_gov.df$all.right==1,]
military_gov.treat.right$num <- rowSums(military_gov.treat.right==1)-2  
prop.table(table(military_gov.treat.right$num))

military_gov.treat.left <- military_gov.df[military_gov.df$all.treat==1&military_gov.df$all.right==0,]
military_gov.treat.left$num <- rowSums(military_gov.treat.left==1)-1  
prop.table(table(military_gov.treat.left$num))

military_gov.control.right <- military_gov.df[military_gov.df$all.treat==0&military_gov.df$all.right==1,]
military_gov.control.right$num <- rowSums(military_gov.control.right==1)-1  
prop.table(table(military_gov.control.right$num))

military_gov.control.left <- military_gov.df[military_gov.df$all.treat==0&military_gov.df$all.right==0,]
military_gov.control.left$num <- rowSums(military_gov.control.left==1)  
military_gov.control.left$num <- factor(military_gov.control.left$num, levels = c(0:4))
prop.table(table(military_gov.control.left$num))
df <- data.frame(waveschosen=c(0:4),percent=NA)

df$percent <- prop.table(table(military_gov.control.left$num))
df$group <- "control left"
df2 <- data.frame(waveschosen=c(0:4),percent=NA)
military_gov.control.right$num <- factor(military_gov.control.right$num, levels = c(0:4))

df2$percent <- prop.table(table(military_gov.control.right$num))
df2$group <- "control right"

df3 <- data.frame(waveschosen=c(0:4),percent=NA)
military_gov.treat.right$num <- factor(military_gov.treat.right$num, levels = c(0:4))

df3$percent <- prop.table(table(military_gov.treat.right$num))
df3$group <- "treat right"

df4 <- data.frame(waveschosen=c(0:4),percent=NA)
military_gov.treat.left$num <- factor(military_gov.treat.left$num, levels = c(0:4))

df4$percent <- prop.table(table(military_gov.treat.left$num))
df4$group <- "treat left"

milgov.df <- rbind (df,df2,df3,df4)

## satisfaction with police

all$pre_inst_police <- ifelse(all$pre_inst_police==1|all$pre_inst_police==0,0,1)
all$inst_police <- ifelse(all$inst_police==1|all$inst_police==0,0,1)
all$inst_police_f1 <- ifelse(all$inst_police_f1==1|all$inst_police_f1==0,0,1)
all$inst_police_f2 <- ifelse(all$inst_police_f2==1|all$inst_police_f2==0,0,1)
all$inst_police_f3 <- ifelse(all$inst_police_f3==1|all$inst_police_f3==0,0,1)


inst_police.df <- data.frame(all$inst_police, all$inst_police_f1,all$inst_police_f2,all$inst_police_f3,all$treat,all$right)
inst_police.treat.right <- inst_police.df[inst_police.df$all.treat==1&inst_police.df$all.right==1,]
inst_police.treat.right$num <- rowSums(inst_police.treat.right==1)-2  
prop.table(table(inst_police.treat.right$num))

inst_police.treat.left <- inst_police.df[inst_police.df$all.treat==1&inst_police.df$all.right==0,]
inst_police.treat.left$num <- rowSums(inst_police.treat.left==1)-1  
prop.table(table(inst_police.treat.left$num))

inst_police.control.right <- inst_police.df[inst_police.df$all.treat==0&inst_police.df$all.right==1,]
inst_police.control.right$num <- rowSums(inst_police.control.right==1)-1  
prop.table(table(inst_police.control.right$num))

inst_police.control.left <- inst_police.df[inst_police.df$all.treat==0&inst_police.df$all.right==0,]
inst_police.control.left$num <- rowSums(inst_police.control.left==1)  
inst_police.control.left$num <- factor(inst_police.control.left$num, levels = c(0:4))
prop.table(table(inst_police.control.left$num))
df <- data.frame(waveschosen=c(0:4),percent=NA)

df$percent <- prop.table(table(inst_police.control.left$num))
df$group <- "control left"
df2 <- data.frame(waveschosen=c(0:4),percent=NA)
inst_police.control.right$num <- factor(inst_police.control.right$num, levels = c(0:4))

df2$percent <- prop.table(table(inst_police.control.right$num))
df2$group <- "control right"

df3 <- data.frame(waveschosen=c(0:4),percent=NA)
inst_police.treat.right$num <- factor(inst_police.treat.right$num, levels = c(0:4))

df3$percent <- prop.table(table(inst_police.treat.right$num))
df3$group <- "treat right"

df4 <- data.frame(waveschosen=c(0:4),percent=NA)
inst_police.treat.left$num <- factor(inst_police.treat.left$num, levels = c(0:4))

df4$percent <- prop.table(table(inst_police.treat.left$num))
df4$group <- "treat left"
police.df <- rbind (df,df2,df3,df4)

p1 <- ggplot(data=pardon.df, aes(x=waveschosen, y=percent)) +
  geom_bar(stat="identity", width=0.5,colour="black") +
  facet_grid(. ~ group) +
  geom_text(data=pardon.df, aes(x = waveschosen, y = (percent + .05), label = paste0(round(percent*100,1),"%")), colour="black", size = 2.5) +
  theme_bw() + theme(axis.text.y=element_blank(),
                     axis.ticks.y=element_blank(),strip.background =element_rect(fill="white")) +
  ggtitle("Support pardoning perpetrators") +
  xlab("Waves chosen") + ylab(NULL)
p2<-ggplot(data=churchtrust.df, aes(x=waveschosen, y=percent)) +
  geom_bar(stat="identity", width=0.5,colour="black") +
  facet_grid(. ~ group) +
  geom_text(data=churchtrust.df, aes(x = waveschosen, y = (percent + .05), label = paste0(round(percent*100,1),"%")), colour="black", size = 2.5) +
  theme_bw() + theme(axis.text.y=element_blank(),
                     axis.ticks.y=element_blank(),strip.background =element_rect(fill="white")) +
  ggtitle("Trust in church") +
  xlab("Waves chosen") + ylab(NULL)
p3<-ggplot(data=govsat.df, aes(x=waveschosen, y=percent)) +
  geom_bar(stat="identity", width=0.5,colour="black") +
  facet_grid(. ~ group) +
  geom_text(data=govsat.df, aes(x = waveschosen, y = (percent + .05), label = paste0(round(percent*100,1),"%")), colour="black", size = 2.5) +
  theme_bw() + theme(axis.text.y=element_blank(),
                     axis.ticks.y=element_blank(),strip.background =element_rect(fill="white")) +
  ggtitle("Satisfaction with government") +
  xlab("Waves chosen") + ylab(NULL)
p4<-ggplot(data=democracy.df, aes(x=waveschosen, y=percent)) +
  geom_bar(stat="identity", width=0.5,colour="black") +
  facet_grid(. ~ group) +
  geom_text(data=democracy.df, aes(x = waveschosen, y = (percent + .05), label = paste0(round(percent*100,1),"%")), colour="black", size = 2.5) +
  theme_bw() + theme(axis.text.y=element_blank(),
                     axis.ticks.y=element_blank(),strip.background =element_rect(fill="white")) +
  ggtitle("Satisfaction with Democracy") +
  xlab("Waves chosen") + ylab(NULL)
p5<-ggplot(data=milgov.df, aes(x=waveschosen, y=percent)) +
  geom_bar(stat="identity", width=0.5,colour="black") +
  facet_grid(. ~ group) +
  geom_text(data=milgov.df, aes(x = waveschosen, y = (percent + .05), label = paste0(round(percent*100,1),"%")), colour="black", size = 2.5) +
  theme_bw() + theme(axis.text.y=element_blank(),
                     axis.ticks.y=element_blank(),strip.background =element_rect(fill="white")) +
  ggtitle("Support for military government") +
  xlab("Waves chosen") + ylab(NULL)
p6<-ggplot(data=police.df, aes(x=waveschosen, y=percent)) +
  geom_bar(stat="identity", width=0.5,colour="black") +
  facet_grid(. ~ group) +
  geom_text(data=police.df, aes(x = waveschosen, y = (percent + .05), label = paste0(round(percent*100,1),"%")), colour="black", size = 2.5) +
  theme_bw() + theme(axis.text.y=element_blank(),
                     axis.ticks.y=element_blank(),strip.background =element_rect(fill="white")) +
  ggtitle("Satisfaction with Police") +
  xlab("Waves chosen") + ylab(NULL)

grid.arrange(p1, p2,p3,p4,p5,p6, nrow = 3)
#g <- arrangeGrob(p1, p2, p3,p4,p5,p6, nrow=3) #generates g


##############################################################
######Table A20. Political insitutions - dropping missing observations
##############################################################

load(file = "all.Rdata")
est.ate<-function(dv, predv, df){
  summary(fit.1 <- lm(dv~df$treat + df$pre_ideology_1 +
                        + predv*df$date_diff + df$base_gender +df$age + df$v))
  vcv <- vcovHC(fit.1)
  n <- nobs(fit.1)
  result <- coeftest(fit.1, vcv)[2, 1:4] / itt.d
  result <- list("point.estimate" = result[1],"se"=result[2],"pvalue"=result[4], "obs" = n)
  return(result)
}

# This estimates ATE when we don't have a pre-treatment measurement
est.ate.np<-function(dv, df){
  summary(fit.1 <- lm(dv~df$treat + df$pre_ideology_1 + df$base_gender +df$age+df$v))
  vcv <- vcovHC(fit.1)
  n <- nobs(fit.1)
  result <- coeftest(fit.1, vcv)[2, 1:4] / itt.d
  result <- list("point.estimate" = result[1],"se"=result[2],"pvalue"=result[4], "obs" = n)
  return(result)
}

# pol inst DVs
dem <- est.ate(all$democracy, all$pre_democracy, all)
mil <- est.ate(all$military_gov, all$pre_military_gov, all)
gov_sat <- est.ate(all$inst_gov, all$pre_inst_gov, all)
mil_sat <- est.ate(all$inst_mil, all$pre_inst_mil, all)
pol_sat <- est.ate(all$inst_police, all$pre_inst_police, all)
gov_trust <- est.ate(all$conf_gov, all$pre_conf_gov, all)
mil_trust <- est.ate(all$conf_mil, all$pre_conf_mil, all)
pol_trust <- est.ate(all$conf_police, all$pre_conf_police, all)
church_trust <- est.ate(all$conf_church, all$pre_conf_church, all)

##############################################################
######Table A21. Transitional justice - dropping missing observations
##############################################################

# TJ DVs
advance <- est.ate.np(all$justice_advance, all)
justice_account <- est.ate.np(all$justice_account, all)
compensation <- est.ate(all$current_recomp, all$pre_current_recomp, all)
judicial <- est.ate.np(all$tj_judicial, all)
inst_apology <- est.ate.np(all$tj_apology, all)
apologize <- est.ate.np(all$policies_apologize, all)
compensate <- est.ate.np(all$policies_compensate, all)
pardoned <- est.ate.np(all$policies_pardon, all)

##############################################################
######Table A22. Balance on measurements collected at baseline among 
# nonparticipants and participants
##############################################################

# Note that not all subjects who eventually participated in our experiment completed the baseline - nonetheless, many did
load(file = "baseline.Rdata")
nonparticipants <- baseline[!(baseline$ID %in% all$ID),]
participants <- baseline[(baseline$ID %in% all$ID),]

t.test(participants$female, nonparticipants$female)
t.test(participants$ideology, nonparticipants$ideology)
t.test(participants$pinochet, nonparticipants$pinochet)
t.test(participants$pinochet_london, nonparticipants$pinochet_london)
t.test(participants$prosecution, nonparticipants$prosecution)
